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lute control of propagated error by ensuring truncation error is 
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modules to control an effective radiation dose at the target of 
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APPARATUS, METHOD AND PROGRAM 
STORAGE DEVICE FOR DETERMINING 
HIGH-ENERGY NEUTRON/ION TRANSPORT 
TO A TARGET OF INTEREST 

5 

Pursuant to 35 U.S.C. §119, the benefit of priority from 
provisional application 60/877,012, with a filing date of Dec. 

11, 2006, is claimed for this non-provisional application. 

ORIGIN OF THE INVENTION to 

The invention described herein was made by employees of 
the United States Government and may be manufactured and 
used by or for the Government for Government purposes 
without payment of any royalties thereon or therefore. 15 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

This invention relates in general to radiation shield designs, 20 
and more particularly to an apparatus, method and program 
storage device for calculating high-energy neutron/ion trans- 
port to a target of interest. 

2. Description of the Related Art 

The capability to make diagnostic assessments of radiation 25 
exposure is needed to support a wide range of radiation expo- 
sure events. Moreover, the question of risk from radiation 
exposure is a much-debated topic of discussion. Every person 
receives daily “background” radiation from a variety of natu- 
ral sources: from cosmic rays and radioactive materials in the 30 
Earth, from naturally occurring radionuclides in food, and 
from inhaling particulate decay products of radon gas. One 
area of increased radiation exposure risk to human results 
from advancing aircraft technology that allows higher oper- 
ating altitudes thereby reducing the protective cover provided 35 
by the Earth’s atmosphere from extraterrestrial radiations. 
This increase in operating altitudes is taken to a limit by 
human operations in space. Space radiation is likely to be the 
ultimate limiting factor for future human deep space explo- 
ration. Understanding the space radiation environment is 40 
essential for risk assessment of orbit/crew selection and pro- 
vides the scientific basis of countermeasures for shielding 
materials (affecting flight weight/cost), radio-protectants, 
and pharmaceuticals. Every tissue/material/part installed on a 
space mission requires radiation risk analysis. While the 45 
present invention is described here with reference to space- 
craft, those skilled in the art will recognize that the principles 
discussed herein and the embodiments of the invention 
described herein are also applicable to other applications and 
industries, such as aircraft design, material development, and 50 
proton cancer therapy. 

The propagation of galactic ions through extended matter 
and determination of the origin of these ions has been the 
subject of many studies. For example, a one-dimensional 
equilibrium solution was proposed early to show that the light 55 
ions have their origin in the breakup of heavy particles. How- 
ever, the one dimensional equilibrium solution did not 
include ionization energy loss and radioactive decay. Later, 
the one-dimensional propagation was shown to be simplistic 
and that leakage at the galactic boundary must be taken into 60 
account. The leakage was found to be approximated as a 
superposition of nonequilibrium one-dimensional solutions. 

A solution to the steady-state equations was given as a Volt- 
erra equation, which was solved to the first order in the frag- 
mentation cross sections by ignoring energy loss. This pro- 65 
vided an approximation of the first-order solution that 
included ionization energy loss and was only valid at relativ- 
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istic energies. An overview of the cosmic ray propagation was 
later provided. A derivation of the Volterra equation included 
the ionization energy loss, but evaluated only the unperturbed 
term. 

These studies focused on only achieving first-order solu- 
tions in the fragmentation cross sections where path lengths 
in the interstellar space are approximately 3 to 4 g/cm 2 . How- 
ever, higher order terms cannot be ignored in accelerator or 
space shielding transport problems. In addition to this sim- 
plification, previous cosmic ray models have neglected the 
complicated three-dimensional nature of the fragmentation 
process. 

Several approaches to the solution of high-energy heavy 
ion propagation that include ionization energy loss have been 
developed during the last 20 years. However, most have 
assumed the straight-ahead approximation and velocity-con- 
serving fragmentation interactions, whereas only a few have 
incorporated energy -dependent cross sections. An approach 
examining a primary ion beam represented the first-genera- 
tion secondary fragments as a quadrature over the collision 
density of the primary beam. An energy multigroup method 
was used in which an energy-independent fragmentation 
transport approximation was applied within each energy 
group after which the energy group boundaries were moved 
according to continuous slowing-down theory. The energy- 
independent fragment transport equation was solved with 
primary collision density as a source and neglected higher 
order fragmentation. The primary source term extended only 
to the primary ion range from the boundary and the energy- 
independent transport solution was modified to account for 
the finite range of the secondary fragment ions. 

An expression was derived for the ion transport problem to 
the first-order (i.e., first-collision) term and gave an analytical 
solution for the depth-dose relationship. The more common 
approximations used in solving the heavy ion transport prob- 
lem were further examined. The effect of conservation of 
velocity on fragmentation and on the straight-ahead approxi- 
mation was found to be negligible for cosmic ray applica- 
tions. Solution methods for representation of the energy - 
dependent nuclear cross sections were derived. The energy 
loss term and the ion spectra were approximated by simple 
forms for which energy derivatives were evaluated explicitly. 
The resulting ordinary differential equations in terms of posi- 
tion were solved analytically. This approximation results in 
the decoupling of motion in space and a change in energy. The 
energy shifts were replaced by an effective attenuation factor. 
Later, the next higher order (i.e., second-collision) term was 
added. The second-collision term was found to be very impor- 
tant in describing 20 Ne beams at 670 A MeV. The three-term 
expansion was modified to include the effect of energy varia- 
tion of the nuclear cross sections. The integral form of the 
transport equation was also used to derive a numerical march- 
ing procedure to solve the cosmic ray transport problem. This 
method accommodated the energy-dependent nuclear cross 
sections within the numerical procedure. Comparison of the 
numerical procedure with an analytical solution of a simpli- 
fied problem validated the solution technique to approxi- 
mately 1 -percent accuracy. Several solution techniques and 
analytical methods have also been developed for testing 
future numerical solutions of the transport equation. More 
recently, an analytical solution for the laboratory ion beam 
transport problem has been derived with a straight-ahead 
approximation, velocity conservation at the interaction site, 
and energy-dependent nuclear cross sections. 

From an overview of these past developments, the appli- 
cations are divided into two categories: a single-ion species 
with a single energy at the boundary and a broad host of 
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elemental types with a broad continuous energy spectrum. 
Techniques, which will represent the spectrum over an array 
of energy values, require vast computer storage and compu- 
tation speed to maintain sufficient energy resolution for the 
laboratory beam problem. In contrast, analytical methods, 
which are applied as a marching procedure have similar 
energy resolution problems. This is a serious limitation 
because a final (i.e., production) high-charge-and-energy 
(HZE) computation method for cosmic ray shielding must be 
thoroughly validated by laboratory experiments. Some 
researchers hope for a single code, which can be validated in 
the laboratory and used in space applications. More recently, 
a Green’ s function has been derived which can be tested in the 
laboratory and used in space radiation protection applica- 
tions. 

Lastly, the problems of free-space radiation transport and 
shielding has been addressed using a high-charge-and-energy 
(HZE) transport computer program, which is referred to as 
the HZETRN program. The HZETRN program (referred to 
herein as 1 995 HZETRN) has been widely used in prior shield 
design verification and validation processes. Additionally, the 
BRYNTRN code, discussed in F. A. Cucinotta, “Extension of 
the BRYNTRN code to monoenergetic light ion beams,” 
NASA TP-3472, 1994, is a baryon transport code used to 
calculate the energy spectrum of secondary nucleons, and has 
been widely used. 1995 HZETRN is described in detail by J. 
W. Wilson et al. in “HZETRN: Description of a Free-Space 
Ion and Nucleon Transport and Shielding Computer Pro- 
gram,” NASA TP-3495, May 1995, which is hereby incorpo- 
rated by reference in its entirety. 1 995 HZETRN is designed 
to provide fast and accurate dosimetric information for the 
design and construction of space modules and devices. The 
program is based on a one-dimensional space-marching for- 
mulation of the Boltzmann transport equation with a straight- 
ahead approximation. The general Boltzmann equation was 
simplified by using standard assumptions to derive the 
straight-ahead equation in the continuous slowing-down 
approximation and by assuming that heavy projectile breakup 
conserves velocity. The effect of the long-range Coulomb 
force and electron interaction was treated as a continuous 
slowing-down process. Atomic (electronic) stopping power 
coefficients with energies above a few A MeV were calcu- 
lated by using Bethe’s theory including Bragg’s rule, Zie- 
gler’ s shell corrections, and effective charge. Nuclear absorp- 
tion cross sections were obtained from fits to quantum 
calculations and total cross sections were obtained with a 
Ramsauer formalism. Nuclear fragmentation cross sections 
were calculated with a semi -empirical abrasion-ablation frag- 
mentation model. An environmental model was also used to 
provide input to the HZE transport computations. 

Nevertheless, improved spacecraft shield design to support 
planned missions to the moon and Mars requires early entry 
of radiation constraints into the design process to maximize 
performance and minimize costs. Of particular importance is 
the need to implement probabilistic models to account for 
design uncertainties in the context of optimal design pro- 
cesses. These requirements need supporting tools with high 
computational efficiency to enable appropriate design meth- 
ods. 

Accordingly, there is a need for an apparatus, method and 
program storage device for calculating high-energy neutron/ 
ion transport to a target of interest. 

It can also be seen that there is a need for an improved 
radiation shield design apparatus, method and program stor- 
age device that implements improvements to the database, 
basic numerical procedures, and algorithms along with new 
methods of verification and validation to capture a well 
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defined algorithm for engineering design processes to be used 
in an early development phase of space exploration shield 
designs. 

5 SUMMARY OF THE INVENTION 

To overcome the limitations described above and to over- 
come other limitations that will become apparent upon read- 
ing and understanding the present specification, the present 
to invention discloses an apparatus, method and program stor- 
age device for determining high-energy neutron/ion transport 
to a target of interest. 

The present invention solves the above-described prob- 
lems by advancing, verifying and validating the transport 
15 codes for calculating high-energy neutron/ion transport to a 
target of interest. The database, basic numerical procedures, 
and computation method are improved. In addition, bench- 
marks are provided for evaluating further problems, for pro- 
viding code portability and for identifying database drift. 

20 A method for calculating high-energy neutron/ion trans- 
port to a target of interest includes: (1 ) defining boundaries for 
a calculation of a high-energy neutron/ion transport to a target 
of interest; (2) calculating the high-energy neutron/ion trans- 
port to the target of interest using numerical procedures 
25 selected to reduce local truncation error by including higher 
order terms and to allow absolute control of propagated error 
by ensuring truncation error is third order in step size, and 
using scaling procedures for flux coupling terms modified to 
improve computed results by adding a scaling factor to terms 
30 describing production of j -particles from collisions of k-par- 
ticles; and (3) providing the calculated high-energy neutron/ 
ion transport to modeling modules to control an effective 
radiation dose at the target of interest. 

In another embodiment of the present invention, a com- 
35 puter program product embodied in a computer readable 
medium and adapted to perform operations for calculating 
high-energy neutron/ion transport across a material of inter- 
est is provided. The operations include: (1) defining bound- 
aries for a calculation of a high-energy neutron/ion transport 
40 to a target of interest; (2) calculating the high-energy neutron/ 
ion transport to the target of interest using numerical proce- 
dures selected to reduce local truncation error by including 
higher order terms and to allow absolute control of propa- 
gated error by ensuring truncation error is third order in step 
45 size, and using scaling procedures for flux coupling terms 
modified to improve computed results by adding a scaling 
factor to terms describing production of j -particles from col- 
lisions of k-particles; and (3) providing the calculated high- 
energy neutron/ion transport to modeling modules to control 
50 an effective radiation dose at the target of interest. 

In a further embodiment of the present invention, a device 
configured to calculate high-energy neutron/ion transport to a 
target of interest is provided. The device includes memory for 
storing data defining boundaries for a calculation of a high- 
55 energy neutron/ion transport to a target of interest; and a 
processor, coupled to the memory, the processor: (1) calcu- 
lating the high-energy neutron/ion transport to the target of 
interest using numerical procedures selected to reduce local 
truncation error by including higher order terms and to allow 
60 absolute control of propagated error by ensuring truncation 
error is third order in step size, and using scaling procedures 
for flux coupling terms modified to improve computed results 
by adding a scaling factor to terms describing production of 
j-particles from collisions of k-particles; and (2) providing 
65 the calculated high-energy neutron/ion transport to modeling 
modules to control an effective radiation dose at the target of 
interest. 
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These and various other advantages and features of novelty 
which characterize the invention are pointed out with particu- 
larity in the claims annexed hereto and form a part hereof 
However, for a better understanding of the invention, its 
advantages, and the objects obtained by its use, reference 
should be made to the drawings which form a further part 
hereof, and to accompanying descriptive matter, in which 
there are illustrated and described specific examples of an 
apparatus in accordance with the invention. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Referring now to the drawings in which like reference 
numbers represent corresponding parts throughout: 

FIG. 1 is a plot illustrating the geometric relations of quan- 
tities relevant to the transport equations derived from the 
coupled linear Boltzmann equations for a closed convex 
domain according to an embodiment of the present invention; 

FIG. 2 is a plot illustrating the range of ions in aluminum; 

FIG. 3 is a plot illustrating the probability of nuclear reac- 
tion as a function of ion type and energy; 

FIG. 4 is a plot illustrating the integral neutron fluence in an 
aluminum shield using the 1995 HZETRN method and the 
present method for a Sep. 29, 1989 solar particle event; 

FIG. 5 is a plot illustrating the integral proton fluence in 
aluminum shield using the 1995 HZETRN method and the 
present method for the Sep. 29, 1989 solar particle event; 

FIG. 6 is a plot illustrating the integral He 4 fluence in 
aluminum shield using the 1995 HZETRN method and the 
present method for the Sep. 29, 1989 solar particle event; 

FIG. 7 is a plot illustrating the integral H 2 fluence versus 
depth in an aluminum shield using the 1995 HZETRN 
method and the present method for the Sep. 29, 1989 solar 
particle event; 

FIG. 8 is a plot illustrating the integral H 3 fluence versus 
depth in an aluminum shield using the 1995 HZETRN 
method and the present method for the Sep. 29, 1989 solar 
particle event; 

FIG. 9 is a plot illustrating the integral He 3 fluence versus 
depth in an aluminum shield using the 1995 HZETRN 
method and the present method for the Sep. 29, 1989 solar 
particle event; 

FIGS. 10 a-d are plots illustrating the total dose and dose 
equivalent for the Webber benchmark SPE spectrum for alu- 
minum and iron on water; 

FIGS. 11 a-b are plots showing numerical errors in proton 
spectra for analytic SPE and GCR benchmarks versus energy 
index; 

FIG. 12 is a plot illustrating a comparison of the results 
derived from the BRYTRN (version 3) method, the 1995 
HZETRN method (including ten years of drift), and the 
present method; 

FIG. 13 is a flow chart of the present method according to 
an embodiment of the present invention; and 

FIG. 14 illustrates a system according to an embodiment of 
the present invention. 

DETAILED DESCRIPTION OF THE INVENTION 

In the following description of the embodiments, reference 
is made to the accompanying drawings that form a part 
hereof, and in which is shown by way of illustration the 
specific embodiments in which the invention may be prac- 
ticed. It is to be understood that other embodiments may be 
utilized as structural changes may be made without departing 
from the scope of the present invention. 
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The present invention provides an apparatus, method and 
program storage device for calculating high-energy neutron/ 
ion transport to a target of interest, and is discussed in J. W. 
Wilson et al. in “Standardized Radiation Shield Design 
5 Method: 2005 HZETRN,” 06ICES-1 8, which is hereby incor- 
porated by reference herein in its entirety. 

Crewmembers in a space module will be exposed to both 
ionizing and non-ionizing radiation. Ionizing radiation, 
which breaks chemical bonds in biological systems, can have 
to immediate (acute) as well as latent effects, depending on the 
magnitude of the radiation dose absorbed, the species of 
ionizing radiation, and the tissue affected. The ionizing radia- 
tion in space is comprised of charged particles, uncharged 
particles, and high-energy electromagnetic radiation. The 
15 particles vary in size from electrons (beta rays) through pro- 
tons (hydrogen nuclei) and helium atoms (alpha particles), to 
the heavier nuclei encountered in cosmic rays, e.g., HZE 
particles (High Z and Energy, where Z is the charge). They 
may have single charges, either positive (protons, p) or nega- 
20 tive (electrons, e); multiple charges (alpha or HZE particles); 
or no charge, such as neutrons. The atomic nuclei of cosmic 
rays, HZE particles, are usually completely stripped of elec- 
trons and thus have a positive charge equal to their atomic 
number. 

25 The ionizing electromagnetic radiation consists of x-rays 
and gamma-rays, which differ from each other in their energy 
and add little to extraterrestrial space exposures. By conven- 
tion, X-rays have a lower energy than the gamma -rays, with 
the dividing line being at about 1 MeV. In general, x-rays are 
30 produced either by the interaction of energetic electrons with 
inner shell electrons of heavier elements or through the brak- 
ing radiation mechanism when deflected by the Coulomb 
field of the atomic nuclei of the target material. Gamma-rays 
are usually products of the de-excitation of excited heavier 
35 elements. 

Mass shielding is the main means of protecting crewmem- 
bers from space radiation. Space modules are constructed 
with an outer skin and associated structural members, and 
sometimes an outer micrometeoroid/space debris shield. In 
40 addition, the space module contains specialized equipment 
with considerable mass and internal structural features (e.g., 
walls, cabinets) which can provide some additional shielding, 
but in only some specific directions as these masses are not 
distributed uniformly and/or isotropically. 

45 Improved spacecraft shield design requires early entry of 
radiation constraints into the design process to maximize 
performance and minimize costs. The atomic and nuclear 
processes associated with space radiation occur over very 
short time scales (microseconds) compared with the secular 
50 variations of the space environment. This allows the use of a 
time independent master equation represented by a steady- 
state Boltzmann description balancing gains and losses of the 
particle fields, e.g., Galactic Cosmic Rays (GCR) and Solar 
Particle Events (SPE), interacting with the shield material 
55 (including the human tissues) . This equation may be reduced 
to a readily soluble numerical process. 

The specification of the interior environment within a 
spacecraft and evaluation of the effects on the astronaut is at 
the heart of the space radiation protection problem. The rel- 
60 evant transport equations are the coupled linear Boltzmann 
equations for a closed convex domain. 

FIG. 1 is a plot 100 illustrating the geometric relations of 
quantities relevant to the transport equations derived from the 
coupled linear Boltzmann equations for a closed convex 
65 domain according to an embodiment of the present invention. 
FIG. 1 establishes the frame of reference for 4>^(x, Q, E) 
representing the flux of ions of type j at x 110 with motion 
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along Q 112, where T 114 is the point on the boundary 
connected to x 110 along Q 112 and n 116 is the unit normal 
vector at the boundary surface at point T 114. The coupled 
linear Boltzmann equations are derived on the basis of con- 
servation principles for the flux density (particles/cm 2 -sr-s- 
A-MeV) <|) .(x Q, E) for particle type j as: 

Q,E)=Z / Jo Jk (Q,Q , £,E , )<bk(x, & ',E')dQ ’dE’- 

Oj(E)§j(x,Q,E), (1) 

where cr^E) and cr^QjQ'jEjE') are the shield media macro- 
scopic cross sections. The a^QjQ'jEjE') represent all those 
processes by which type k particles moving in direction Q' 
with energy E' produce a type j particle in direction Q with 
energy E (including decay processes). Note that there may be 
several reactions that produce a particular product, and the 
appropriate cross sections for equation (1) are the inclusive 
ones. Exclusive processes are functions of the particle fields 
and may be included once the particle fields are known. Note, 
at times Q.(x, Q, E) will be loosely referred to as either flux or 
fluence and the usage should be clear from the context. The 
time scale of the processes in equation (1) are at most on the 
order of microseconds while the time scales of boundary 
conditions are on the order of minutes or longer, leaving the 
resulting interior fields in equilibrium with the particles at the 
boundary. 

The total cross section cr(E) with the medium for each 
particle type is: 

°j(E) = Oj at (E)+o/ l (E)+Oj r (E), (2) 

where the first term refers to collision with atomic electrons, 
the second term is for elastic scattering on the nucleus, and the 
third term describes nuclear reactions where the minor 
nuclear inelastic processes (excited single particle states) 
have been ignored except for low energy neutron collisions. 
The corresponding differential cross sections are similarly 
ordered. Many atomic collisions (~10 6 ) occur in a centimeter 
of ordinary matter, whereas ~10 3 nuclear coulomb elastic 
collisions occur per centimeter, while nuclear scattering and 
reactive collisions are separated by a fraction to many centi- 
meters depending on energy and particle type. The o^E) 
term includes the nuclear decay processes. Solution methods 
first use physical perturbations based on the ordering of the 
cross sections with the frequent atomic interactions as the first 
physical perturbation with special methods used for neutrons 
for which atomic cross sections are zero. The first physical 
perturbation to be treated is the highly directed atomic colli- 
sions with mean free paths on the order of micrometers as 
observed in nuclear emulsion. 

FIG. 2 is a plot 200 illustrating the range of ions in alumi- 
num. The usual approximation is the continuous slowing 
down approximation leading to well-specified range-energy 
relations as shown in FIG. 2. In FIG. 2, the range 210 is 
plotted against the energy 220 of ions in aluminum for a range 
of Z values 230. In FIG. 2, the energy straggling is neglected. 
This energy straggling will be discussed later. The next term 
is the highly directed multiple Coulomb scattering. This term 
is usually neglected in many models, but is of great impor- 
tance in understanding the transport of unidirectional ion 
beams leading to beam divergence. The remaining nuclear 
reactive processes have been given main attention in past 
code developments. 

Continuous Slowing Down Approximation 

The collisions with atomic electrons preserve the identity 
of the ion and the differential cross sections are given as: 

CT y ./'(Q,Q'A£>2„o^%(£')6(Q-Q')6 / -*x6(£+4- 1 e„- 

E'\ ’ (3) 
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where n refers to the atomic/molecular excited states with 
excitation energies € w including the continuum. Note, the 
factor A~ l results from the units of E of A MeV (equivalent 
unit of MeV/nucleon with atomic weight A-). Although the 
5 atomic/molecular cross-sections o at jn (E) are large (~10 -16 
cm 2 ), the energy transfers €„ are small (^ 1-1 00 eV) compared 
to the particle energy. The atomic/molecular terms of equa- 
tion (1) may be written as: 

10 

e, E')<p k (x, el', E')dEi'dE' - o-y (£>£/(*, el, ' 4 1 

E = Z n (Ej,n(E + Aj l s n )(f>j(x, EL, E+ A] l s n ) - of (E)(f> j(x, EL, E) = 
^ n {(Ej, n {E)(f)j{x, EL, E) + Afs n d E [oJ^EWjix, EL, £)]} - 
v"j{E)<t>j{x , EL,E) + 0(e 2 n ) = 

d E [Sj(E)(Pj(x, a £)] + <?(£*), 

20 where the stopping power S^.(E) is given as the sum of energy 
transfers and atomic excitation cross sections as: 

(5) 

25 The higher order terms of equation (4) are neglected in the 

continuous slowing down approximation (csda). Evaluation 
of the stopping power by equation (5) is deceptively simple in 
that all of the excited states including continuum states of the 
atomic/molecular system need to be known. Furthermore, the 
30 projectile remains a bare ion except at low energies, where the 
projectile ion atomic orbital states begin to resonate with the 
electrons of the media leading to electron capture and lower- 
ing of the ion charge. Equation (1) can be written in the csda 
as: 

Q • V<J)/x, Q,E)-Ar l d E [Sj{E)§j(x, Q,£)]=2ja^(Q,Q'^, 

Q \E)dQ 'dE'-Oj(E)<bj(x, QE), (6) 

where the right-hand side of equation (6) excludes the atomic/ 
molecular processes now appearing on the left as an energy 
40 shifting operator in addition to the usual drift term. Neutral 
particles would have null atomic cross sections for which the 
stopping term of equation (6) does not appear. Application of 
csda in both laboratory and space shielding has been wide- 
spread, including the resulting errors. Equation (6) can be 
45 rewritten as an integral equation: 

(E')]Q,Q'F')}/Sj(E)Pj(E), (7) 

where, again referring to FIG. 1, T 114 is the point on the 
50 boundary connected to x 110 along Q 112, E Y =R 7 -1 [p-d+R / ], 
p is the projection of x 110 onto Q 112, d is the projection of 
T 114 onto Q 112, R,(E) is the distance an ion of type j of 
energy E will travel before losing all of its energy to excitation 
of atomic electrons, and P^(E) is the probability a type j ion of 
energy E will have a nuclear reaction in coming to rest in the 
media. The usual range-energy relation is given by: 

Rj(E)=UjdEVSj(E) (8) 

60 FIG. 3 is a plot 300 illustrating the probability of nuclear 
reaction 310 as a function of ion type 320 and energy 330. 
With reference to FIG. 3, the nuclear attenuation function is 
given by: 

65 P J (E)^xp[-Uj° r j(E')dEYS j (E)], (9) 

where the integral domains in equations (8) and (9) extend 
over the lull energy range {0, E}. 
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Straight-Ahead Approximation 

The approach to a practical solution of equation (7) is to 
develop a progression of solutions from the simple to the 
complex, allowing early implementation of high-perfor- 
mance computational procedures and establishing a converg- 5 
ing sequence of approximations with established accuracy 
criteria and means of verification. The lowest order approxi- 
mation using the straight-ahead approximation uses the 
Monte Carlo methods, in which the differential cross sections 
are approximated as: 10 

( 10 ) 

resulting in dose and dose equivalent per unit fluence to be 
within the statistical uncertainty of the Monte Carlo result 
obtained using the fully angle dependent cross sections. The 15 
relation of angular dependent cross sections to spacecraft 
geometry in space application is examined using asymptotic 
expansions about angular divergence parameters demonstrat- 
ing errors in the straight-ahead approximation to be on the 
order of the square of the ratio of distance of divergence to 20 
radius of curvature of the shield (a small error in most space 
systems). 

Equations (6) and (7) were examined for HZE ions using 
the following form for the projectile fragmentation cross 
sections as: 25 

o r Jk (Q,Q%E , )=o r Jk (E')N t exp{- [£WE-QVE] 2 /2e jk 2 }, (1 1) 

where o J C Ar (E t ) is the cross section for producing fragment j 
from ion k, N r is the normalization constant for the exponen- 
tial function, and € Jk i s the momentum dispersion parameter in 30 
the reaction. Substituting the interactive form of equation 
(11) into the integral term of the Boltzmann equation (6) 
yields 

Q \E')dQ tf£'=2o^(£') Q, 

E)+E3^> k (x, Q ,E)V [^l(2mE)]+Q+d^k(x, QE) 35 

e//(2 mE)}, (12) 

where the second term on the right hand side of equation (12) 
results from corrections in assuming the velocity of the ion is 
preserved in the interaction, and the third term is error result- 
ing from the straight-ahead assumption. The surprising result 
is that the velocity conserving assumption is inferior to the 
straight-ahead approximation for the nearly isotropic space 
radiation. Under approximations examined in equations (4) 
and (12), there are great simplifications in the Boltzmann 
equation, as given below 

Q • V (^(x, Q,E)-Aj- l d E [Sj(E)tyj(x, Q,K)]=L&‘ J ^E)^ k (x, 

Q ,E)-<fj{E)tyj(x,&,E) (13) 

which is strictly applicable to the HZE ions (Z>2). The light 
ions and neutrons have additional complications arising from 50 
the broad energy spectra associated with their production, 
although the more favorable straight-ahead approximation is 
useful, as indicated in equation (12). The corresponding light 
ion (and neutron) Boltzmann equation is: 

55 

(x, Q ',E')dE'-Gj(E)§j{x, Q,E), (14) 

where the straight-ahead approximation as given by equation 
(10) is used. Equations (13) and (14) have sufficient simplic- 
ity to allow an approach for both space and laboratory appli- 60 
cations. The main force of the laboratory applications allow 
detailed model testing of the many atomic/molecular and 
nuclear processes. 

Marching Procedures and HZETRN 

Both the 1 995 HZETRN and 2005 HZETRN are based on 65 
the solution of equation (14) using straight-ahead approxima- 
tion, as described by equations (10) through (12). Specializ- 
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ing the solution along a ray Q in the direction of the x-axis 
results in the one dimensional description of the Boltzmann 
equation as: 

[dx-Aj-^S^+o^/x^Hj Gj k (E,E)<k k (x,E)dE ' (15) 

Where a Jk (E^E) are approximated for nucleons. An immedi- 
ate problem is the near singular nature of the differential 
operator, and transformation from energy to residual range 
coordinates as used in developing the Green’s function 
greatly relieves this problem. Unlike the Green’s function 
development, numerical procedures are simplified by intro- 
ducing only a single residual range coordinate for all ions. 
The residual proton range r is used as the common coordinate: 

r=$ 0 E dE'/S(E) (16) 

and the residual range of other particle types is related 
through a scaling parameter y J =Z J 2 IA J . as r, wherein A . 
and Zj are mass number and charge number, respectively, 
which fails at low energies corresponding to low residual 
range due to electron capture into atomic orbitals character- 
istic to each ion type. The corresponding transport equation 
is: 

[3*- v A +0>W]^A^) = 2j’ r *(vyv^ Jfc (i; r')y k (x,r')dr\ (17) 

where scaled flux is now for neutral particles such as 
neutrons are taken as unity in scaling relations): 

^j{x,r)=VjS(E)%(x,E), (18) 

and the sealed differential cross sections are: 

Sj k {r,r')=S{E)v jk {E,E). (19) 

Errors in scaling of proton-stopping and range parameters 
in arriving at the approximate transport equation (17) are 
compensated in part by solutions of equation (17) approach- 
ing a low energy equilibrium spectrum for ions given by: 

VjS{E)tyj(x,E) ^constant, (20) 

where the constant is fixed by the higher ion energy. In dis- 
tinction, the solution to equation (15) for ions has the low 
energy equilibrium spectrum: 

Af 1 Sj(E)ty J (x,E) ^constant, (21) 

which is also fixed by the higher energy flux for which the 
range scaling relation vyr^r has better validity and the two 
constants are nearly equal so that equation (21) has improved 
accuracy over equation (20) at lower energies. This fact 
requires alteration of the flux unsealing relations as 
demanded by equation (21) to maintain accuracy at the lower 
energies. From equations (20) and (21), the simplicity of 
numerically solving equation (17) can be understood over a 
numerical solution based on equation (15). The solution to 
equation (17) approaches a constant at small residual ranges, 
allowing large separations in r grid values with smooth 
extrapolation to zero range, whereas solutions of equation 
(1 5) vary as the nearly singular l/S^E) for which small E grid 
spacing is required, leading to slow computational proce- 
dures. The assumptions in equation (17) are tested and 
unsealed according to relation (21) as shown later herein. 

The confusion caused by different scaling methods and 
associated coordinates for numerical procedures is justified 
by the simplification of the numerical representation of flu- 
ence of all particle types over a common residual range grid 
and simplification of the numerical procedures leading to 
high performance codes. Still a straightforward finite differ- 
encing of equation (17) can introduce unstable roots, as had 
plagued the thermal transport problem for many years. The 
differential operator of equation (7) is inverted as shown by: 

r )=exp [-y sxflip/o, r+VjX)+Lj 0 x $ r+vjx E ex p 

*')] ( v / v k) s jk( r +v jx ' r')xty k (x-x \ r')dr'dx \ (22) 
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where the exponential is the integrating factor related to 
attenuation of the j type ions with: 

^(r,x)=$ x Oj(r+vjx')dx', (23) 

which is related to equation (9). Equation (22) is a Volterra 5 
equation and can be solved either as a Neumann series or with 
marching procedures. Note that the inverse mapping is taken 
as: 


$j{x,E)=Aj-tyj(x,r)iSj{E), (24) 

10 

to guarantee the equilibrium solution given as equation ( 21 ) at 
low energies away from the boundaries (note, the proton 
stopping power is used in case of unsealing the neutron flux). 
The equilibrium constant resulting from equation (22), and 
given in equation ( 20 ), is assumed to differ little from condi- 15 
tion (21), for which the inverse mapping of equation (24) is 
most accurate. These approximations are verified later herein. 

Two tracks are taken in implementing a marching proce- 
dure for equation ( 22 ) depending on particle type as 
demanded by the character of the nuclear processes. The 2 q 
problem naturally divides into “light ions,” which will refer to 
all ions with atomic mass of four or less including neutrons, 
and into high charge-energy (HZE) ions having atomic mass 
greater than 4. The distinction arises from the energy and 
angle distributions of the double differential cross sections, 
for which the HZE ions leaving a projectile fragmentation 
event have velocity nearly equal to that of the projectile, as 
approximated by equation (11). Although the light ions are 
assumed to travel in the same direction as the projectile (see 
equation 10 ), they cover a broad energy distribution that 30 
cannot be ignored. The marching procedure is obtained by 
first considering equation ( 22 ) evaluated at x+h, where h is the 
step size as follows: 

% (x+h, r)=exp [-l(r,h ) (x, r+vfa+Hj 0 h f r+vJx ~ exp [~lj 

(r,x')]{v/v^)Sj k (r+vjK ' r')xty k (x+h-x ' r')dr 'dx (25) 35 

Equation (25) may be used to develop a marching step from x 
to x+h once a means to approximate the field function q^(x,r) 
across the subinterval {x, x+h} is provided. If his sufficiently 
small such that 

40 

Ofj(r)h« l, (26) 

then, following lowest order perturbation theory: 

x[)^+^-x;r , )=expE^(r'A-x0]ip*A,r'+v Ar (A-x , )]+O 

(h), (27) 

45 

which may be used to approximate the integral in equation 
(25), giving results for the fields 0(h 2 ) as required to control 
the propagated error. Substituting equation (27) into (25) and 
evaluating the attenuation factors at the interval midpoint 
(mean value theorem) results in: 

ijj j (x+h, r)=exp [-£,(?; Zz) 1 y j(x, r +Vjh ) +2*exp [-C,j(r,h/1)- 

X>k( r> ’ h/2 ) ] x/ r+ v/a/2 00 ( v / v k)F A jk(h, K r '+Vkh/ 

2 )]dr'+0(h 2 ), 

where the integrand has been simplified using 


Ff k (h,r,r')= s Jk (r + Vjx',r'dx' =Fj k (r + Vjh,/)-F jk {r,r'), 

Jo 

and 

ns(r) 

Fjkir,F)= <r jk (E",E')dE", 

Jo 

with e(r) being the energy associated with proton residual 65 
range r, and E-€(r'). Note that if j corresponds to a neutral 
particle, such as the neutron (j=n), then the above expressions 


50 

(28) 

55 

(29) 

60 

(30) 


are evaluated in the limit as v 7 - approaches zero in the range 
scaling relations, resulting in the following (whereas the flux 
scaling factor for neutrons assumes v„=l): 


T|) n (x+h, r)=exp [-ct„ (r)h ]ty n (x, r)+Z k ^„exp [-o n (r)h/2- 
lk(r', Zz/2)]ZzxJ> ( l/v k )s nk (r, r')^ k (x, r'+v k h/2)dr'+ 
exp [- a„ (r)h/2 -o„(r')h/2 ]h$™s nn (r,r')ty n (x, t*)dr', (31) 

and similarly for the neutral k term (k=n) when the j particle 
is charged: 


%(x+h, r)=exp [-’C j (r,h)]^ j (x, H-v/O+S^exp [~lj(r,h/ 

2)-^*(^( h/2)y™(vyV]^^ Jk (h, r, r '+Vjh/2)'xty k [x, r 1 + 
(v J +-v^h/2)]dr'+^p[-X s (r,h/2)-oJ/)h/2]x 

}™v^ A jJJi,r,r'+Vjh/2)\\)Jx,r'+Vjh/2)dr', (32) 

where v w in the flux scaling relation (24) is taken as unity. 
Equations (31) and (32) are solved on an equally space x-grid 
Ax=h apart and a logarithmic spaced r-grid on two subinter- 
vals. The remaining integrals in these equations are approxi- 
mated by: 


r')dr', (33) 

where 00 denotes a chosen upper limit tailored to the specific 
boundary condition. Note that the matrix of K-values can be 
evaluated once on the r-grid and stored for subsequent steps, 
providing high computational efficiency. Equations (31) and 
(32) provide the basis of the light ion transport of both the 
HZETRN 1995 and the BRYNTRN codes. The HZE ion 
projectile (A>4) coupling to the light fragments is contained 
in equations (28) to (32). 

The HZE fragments are produced with nearly the same 
velocity as the projectile ion, as expressed in equation (13), 
and results in the simplified Boltzmann equation: 

for which the scaled equations result in contributions from all 
HZE ions (with A^>4) as: 


r ) = exp [-^x)]i|),.(0, r+vyO+Sj ^exp K//A)] (v/ 
v*)xcr^(r+v/)%(x-x ( r+vjx')dx'. (35) 

The corresponding marching equation is given as: 

il>/*+/z,r)=exp [-lj(r,h)]^)j(x, r+v/)+2j 0 ; 'exp[-^.(^x 1 )] 

(Vj/v k )xOj k (r+VjX')^ k (x+h -x \ r+Vjx)dx \ (3 6) 

for which the integrand can be approximated for sufficiently 
small h using: 

Tj^x+Zz-x r+vA) = exp [-Z, k (r+VjX '/z-x')]xt|)^/x, r+vjx ’+ 

v k (h-x')]+0(h), (37) 

allowing the following simplification: 

y,-(x+/z,r)=exp [-^(nh)]^, r+vJi)+J.j$ 0 h Qxy[-X ) j(r,x')- 
Xskir+vjx \ h-x')] (v/v k )x 0 j k (r+vjx')^ k [x, r+vjx '+v k 
(h-x')]dx'. (38) 

To evaluate equation (38), the mean value theorem that 
guarantees linear terms of the final integral to be zero is used. 
First, the attenuation factor is expanded as: 


Xj(r,x')=$/Oj(r+v/')dx"~$/[Oj(r+Vjhn)+d r Oj(r+Vjh/ 2) 

Vj(x"-h/2)]dx", (39) 

and similarly for: 


l k (r+v J x , ,h-x')=f 0 h x 'o k [r+vjx'+v k (h-x")]dx"~$ 0 h x '[Oj 
(r+VjX '+v k h/2)+d r o j (r+vjx '+v k h/l)v k (x "-h/2)] dx", 


( 40 ) 
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while applying the mean value theorem to the remaining 
factors of equation (38) and neglecting all but linear expan- 
sion terms in the integrand yields: 


ipj{x + h,r) = exp [-£(/■, hftiffj(x, r + vjh) + 

Zk ( vj / v k )o-jk(r+ vjh 1 2)ip k [x,r + (vj + v k )h/ 2] X 

J' exp {-ctj (r + vjh 1 2)x ' - 

o- k [r + (vj + v k )h/2)(h-x')]}dx' = 
exp [-£/(r, h)]ipj{x, r+Vjh) + Z k (vj/v k )o- jk (r + 
vjhl2)ip k [x,r + 

( vj + v k )h/ 2] x [exp{-cr/r + Vjh/2)h} - 
exp{-o"jt [r 4- (vj + v k )h I2)]h}] j {<r k [r + 

(vj + v k )hf2)\ - cr j{r + vjh/2)} + 0(h 2 ), 


to be compared with the 1995 HZETRN algorithm to 0[(v 7 - 
v A )h] given as: 

^j(x+h,r)~exp v/)+ 2 *(v/^)a,*(r)x% 25 

(x, r+v/i){exp [-a/r)/z]-exp [-o k (r)h]}/[o k (r)-Oj 
(r) ]• (42) 

In earlier versions of BRYNTRN for proton/neutron trans- 
port, the flux scaling relation was taken correctly as: 

30 

^(x,r)=S(E)^(x,E), (43) 

but carried over to the latest BRYNTRN for light-ions/neu- 
tron transport. In coupling to HZETRN with scaling given by: 

(44) 35 

there is an inconsistency in flux scaling which must be 
accounted. The appropriate coupling is given in equations 
(38) through (42) with the added factor of v/v* in the field 
coupling terms. The main effects on solution of the Boltz- 
mann equation are expected for the light ions of H 2 , H 3 , and 40 
He 3 with only minor effects on the major light-ion/neutron 
components (n, H 1 , He 4 ). To evaluate these differences in flux 
scaling, the algorithm of equations (31) through (33) have 
been used for comparison with the original light-ion/neutron 
propagator. A 29 Sep. 1989 solar particle event spectrum is 45 
used because of its relation to the 23 Feb. 1956 event repre- 
sented by the proton spectrum (p/cm 2 -MeV) at the boundary 
approximated above 30 MeV by: 

(j>A0,£)=(2-034xl0 7 /(3)x[p(£)/ j p(30)]- 4 - 5 , (45) 50 

where p(E) is the proton momentum (MV) given as: 

p(E)=V[E(E+ 1876)], (46) 

and p is the proton speed relative to the speed of light. A low 55 
energy correction below 30 MeV mainly affecting transport 
results for depths less that 1 g/cm 2 in most materials is also 
added as: 

^(0,£)=1.416xl0 8 xexp[-^(£)/102.1 18]x(£+938)/^ 60 

(E), (47) 

which is in agreement with spectrometer data of the GOES 
satellite. 

FIG. 4 is a plot 400 illustrating the integral neutron fluence 
in an aluminum shield using the 1995 HZETRN computation 65 
method and the present method due to a Sep. 29, 1989 solar 
particle event. In FIG. 4 , the integral fluence 410 , in particles/ 


cm 2 , is plotted against the depths 420 , i.e., g/cm 2 , for both the 
1995 HZETRN computation method 440 and the present 
method 442 . 

FIG. 5 is a plot 500 illustrating the integral proton fluence 
in aluminum shield using the 1995 HZETRN computation 
method and the present method for the Sep. 29, 1989 solar 
particle event. The integral fluence 510 , in particles/cm 2 , is 
plotted against the depths 520 , i.e., g/cm 2 , for both the 1995 
HZETRN computation method 540 and the present method 
542 . 

FIG. 6 is a plot 600 illustrating the integral He 4 fluence in 
aluminum shield using the 1995 HZETRN computation 
method and the present method for Sep. 29, 1989 solar par- 
ticle event. Again, the integral fluence 610 , in particles/cm 2 , is 
plotted against the depths 620 , i.e., g/cm 2 , for both the 1995 
HZETRN computation method 640 and the present method 
642 . 

In FIGS. 4 - 6 , the integral fluence values above 0.01 A MeV 
for neutrons, H 1 , and He 4 with \j= 1 are nearly unchanged, and 
are indistinguishable in FIGS. 4 - 6 , as they are the major 
components produced in reactions and H 1 is dominated by the 
fluence at the boundary over the first half of the mean free 
path. 

FIG. 7 is a plot 700 illustrating the integral H 2 fluence 710 
versus depth 720 in an aluminum shield using the 1995 
HZETRN method 740 and the present method 742 based on 
the Sep. 29, 1989 solar particle event. FIG. 8 is a plot 800 
illustrating the integral H 3 fluence 810 versus depth 820 in an 
aluminum shield using the 1995 HZETRN method 840 and 
the present method 842 based on the Sep. 29, 1989 solar 
particle event. As can be seen in FIGS. 7 - 8 , the H 2 and H 3 
integral fluences are decreased according to their v ; . factors 
with values of Vi and A respectively. 

FIG. 9 is a plot 900 illustrating the integral He 3 fluence 910 
versus depth 920 in an aluminum shield using the 1995 
HZETRN method 940 and the present method 942 based on 
the Sep. 29, 1989 solar particle event. As can be seen in FIG. 
9 , the He 3 integral fluence 910 is increased by the factor of 
v / =56. It is expected that dose will change little as the excess 
of doubly charged He 3 contribution will largely cancel the 
singly charged H 2 and H 3 deficit contributions (approxi- 
mately by a factor of ( 4 A-%) times the total minor contribu- 
tor’s dose). 

The second correction to the propagator algorithm derived 
above, concerns the added accuracy of the HZE propagator to 
0(h 2 ) in equation (41) as opposed to the 1995 HZETRN with 
error term Offy-v^h]. The improved HZE propagator of 
0 (h 2 ) allows control of the propagated error as well as reduc- 
ing the local truncation error as will be demonstrated below. 
Numerical Analysis of Marching Procedures 

There are two variables for which numerical approxima- 
tion enter into the propagator algorithms. The first is in the 
position variable x and the second is the residual range vari- 
able r. The coupling integrals of the Boltzmann equation 
involve integrals over energy that become principally inte- 
grals over residual range for the scaled flux equations, 
although the energy shift operator of the Boltzmann equation 
couples residual range shift and position drift operators along 
the characteristic curves of the transport solution. The prin- 
cipal concern is the necessary control of local truncation 
errors to insure that propagated error is controlled. In consid- 
eration of how errors are propagated, the error introduced 
locally by evaluation of \f(x, r z +h) over the range (energy) 
grid with which it is defined is: 

t|)(x+/z, rf) =exp (- ct/z )ip (x, r t +h ) , (48) 
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whereas the local truncation error is given by: 

t|> (x, r{*-h)=ty int (x, r f +/z)+e f (/z, r t ) . (49) 

After the k th step from the boundary, the numerical solution 
is 

(kh, ^-) =exp (- oh)ty int [(k- 1 )h, r f +/z7+2 x=0 * _1 exp [-a 

(50) 

If the local truncation error is bounded above such that 
€ > (h,r z .)<€(h) for all X, then the propagated error is bounded 
by: 

e pr i ?(^) = 2x=o* _1 exp [- a(£- X)/7 e K (h)<e (h )H K=0 k ~ 1 exp 

[-o(k-X)h 7 «e(/z) [1 -exp (~okh)]/ho, (5 1) 

which is well behaved for all k and h if the local truncation 
error is bounded above by at least 0(h 2 ). The propagated error 
grows to a maximum of €(h)/ha requiring the 0(h 2 ) limitation 
on the local error. The asymptotic bound for deep penetration 
is found to be: 

e prp (h)<e(k)Qxp(-oh)/[l-exp(-oh)], (52) 

emphasizing again the need to control the local truncation 
error as hcr=>0. Earlier BRYNTRN and HZETRN propagator 
algorithms marginally met these requirements. In the reduc- 
tions leading to equations (31), (32) and (41), the error terms 
are 0(h 2 ) when the base algorithms are obtained, but the 
errors associated with the numerical approximation of the 
remaining functions of residual range (or energy) have been 
left so-far unspecified and were the subjects of prior studies. 

Earlier methods assumed approximate log-linear depen- 
dence of all discretized field functions of residual range that 
are on 0(A 2 ) for galactic cosmic ray like spectra, where A is 
the order of the residual range spacing but only O(A) for most 
model solar particle events or trapped proton spectra. 

The original range-grid was derived using a uniform log 
(E)-grid of thirty points converted to range using range-en- 
ergy relations of the transport media. A previous study used a 
90-point log(E)-grid as standard for evaluation of errors in the 
original 30-point grid and a 60-point grid. Maximum errors 
were first quantified to be a few percent in dose and dose 
equivalent at the largest depths of 1 50 g/cm 2 in air. A system- 
atic study of grid generation and numerical interpolation was 
completed. It was found that a uniform log(r)-grid of 
60-points gave an accurate interpolation (fraction of a percent 
of flux) with a fourth order Lagrange interpolation. It was 
desirable at that time to minimize the number of grid points as 
computational time is dominated by evaluation of the integral 
coupling terms and increases as N 2 . It was clear that only the 
midrange errors were significant, so the fully uniform grid 
was replaced with a uniform grid over two sub-domains, 
allowing even greater accuracy with only 30 grid points. An 
excess number of points over the range of 1 g/cm 2 , with fewer 
points at the lower range values, is sufficient. The errors due 
to the residual range grid below 1 g/cm 2 have no effect on the 
propagated error as the step size is on the order of 1 g/cm 2 so 
that this low energy part of the spectrum is deposited in the 
sub -range of the next step. This is facilitated by the scaled flux 
that approaches a constant at these lower energies [see equa- 
tions (20) and (21)]. 

Aside from the issue of numerical interpolation and direct 
effects on the propagation routines, the evaluation of integrals 
of field quantities relates to coupling terms. Past methods 
used the assumed log-linear dependence and evaluated quan- 
tities analytically, arriving at computationally efficient pro- 
cedures (an important feature on contemporary machines at 
that time). Studies of numerical integration errors were made 
using the 90-point solutions as a standard for which the origi- 
nal algorithms for integral flux resulted in errors of less than 
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0.5 percent. It was found that substitution of a three-point 
Simpson’s rule reduced the integration errors by approxi- 
mately an order of magnitude using midpoint values of the 
improved interpolation algorithm with the modified uniform 
5 log(r)-grid on two sub-domains. The reformulated propaga- 
tion routines were found to have a fraction of percent error 
over the transport domain to 1 50 g/cm 2 depths. In every case 
so far studied, the approximations in equation (41) are 
assumed correct and attention is given to evaluation of the 
10 right hand side without reference to the original integral on 
the left side of equation (32). 

The step size convergence within the BRYNTRN algo- 
rithm was examined using the aforementioned modifications 
15 with the 30-point converged results. The step size was varied 
from 1 g/cm 2 to 0. 1 g/cm 2 for which dose for protons con- 
verged quickly but neutrons more slowly. The compromise 
step of 0.5 g/cm 2 is now standard in the BRYNTRN code and 
in the light ion propagator of HZETRN. The current version, 
20 so configured as discussed above with 30 log(r)-grid points, 
results in 5 percent accuracy to 1 50 g/cm 2 and is sufficient for 
most applications. Even so, standard practice now uses 80 
such grid points assuring even improved accuracy for both 
GCR and SPE applications. Furthermore, the number of grid 
25 points is further adjusted to accommodate the simulation of 
geomagnetic cutoff effects while maintaining high numerical 
accuracy. 

Evaluations were made of dose and dose equivalent (as 
given by both the International Commission on Radiological 
30 Protection ICRP 26 and ICRP 60 quality factors) in 30 cm of 
water behind a 20 g/cm 2 shield of aluminum (and alternately 
iron) for the approximation of the 23 Feb. 1956 spectrum 
(p/cm 2 -MeV) given as a P o =100 MV spectrum with 10 9 
protons/cm 2 above 30 MeV in the following: 

35 

^(0, £)=10 9 xexp{ [239. l-p(E)]/l 00 }x(2£+ 1876)/ 

[200xp(£)], (53) 

and comparing with the Monte Carlo results and more mod- 
em Monte Carlo codes using ICRP 60 quality factors. The 
40 present method was evaluated with the ICRP 26 quality fac- 
tors. FIGS. 10 a-d are plots 1002 , 1004 , 1006, 1008 illustrat- 
ing the total dose 1010 and dose equivalent 1020 (ICRP 26) 
for the Webber benchmark SPE spectrum for aluminum 
(FIGS. 10 a-b) and iron (FIGS. 10 c-d) on water. 

45 Testing has been performed with a benchmark by neglect- 
ing the integral term of equation (32) and boundary condition 
given by equation (53) in both the analytical solution and 
1 995 HZETRN code. The analytical solution is given in equa- 
tion (35), neglecting the integral term and unsealing the result 
50 according to equation (24). The initial testing of the present 
method chosen at random from various copies revealed that 
the light-ion/neutron cross section routines were corrupted. 
These were replaced by more accurate (and uncorrupted) 
routines. Now, the transported flux is generally within 1 per- 
55 cent of the analytic solution as is the dose using Simpson’s 
rule, but dose equivalent was found to be low by a few percent. 
Replacing Simpson’s mle by a ten-point Gauss-Legendre 
quadrature brings dose equivalent to within 0.15 percent of 
the analytic result and Gauss-Legendre quadrature will be a 
60 permanent feature of the revised HZETRN computation 
method with comparisons in Table 1 . 

Table 1 shows the comparison of dose and dose equivalent 
(ICRP 60) of penetrating protons from analytical solution and 
the numerical solution (in parenthesis), the comparison of 
65 dose and dose equivalent is shown in Table 1 at various depths 
in water for the analytic benchmark of a Webber spectmm on 
20 g/cm 2 of iron shielding 30 cm water. 
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TABLE 1 


Depth, cm 

Dose, cGy 

Dose equivalent, cSv 

0 

8.405 (8.405) 

11.520 (11.505) 

5 

4.083 (4.074) 

5.009 (4.979) 

10 

2.321 (2.316) 

2.817 (2.800) 

15 

1.417 (1.414) 

1.707 (1.696) 

20 

0.909 (0.907) 

1.089 (1.082) 

25 

0.604 (0.603) 

0.720 (0.716) 

30 

0.412 (0.411) 

0.490 (0.487) 


FIGS, lla-b are plots 1102 and 1104 showing numerical 
errors 1110 and 1112 in proton spectra for analytic SPE (FIG. 

11 a) and GCR (FIG. 116 ) benchmarks versus energy index 
1020 and 1122 . Indexed energies for SPE range from 0.01 to 15 
900 MeV. Indexed energies for GCR range from 0.01 to 
50,000 MeV. From the plots of FIGS. 11 a-b, the percent 
differences 1110 and 1112 of the analytical proton flux and 
the numerically generated proton flux at the iron-shield/water 
interface 1140 and 1142 and at exit of the water slab 1150 and 20 
1152 may be determined. 

The results derived from the plots of FIGS, lla-b provide 
a direct test of the basic propagator methodology, and show 
that the basic propagator methodology is quite accurate. In 
addition to allowing evaluation of the accuracy of basic trans- 25 
port procedures and the nuclear attenuation factors, this 
benchmark provides a direct test of using equation (24) for 
unsealing the numerical solution developed on scaling rela- 
tion (44) and demonstrating the requirements for the low 
energy equilibrium solution of equation (15) to be accurately 
maintained by the approximate numerical propagation 
method. The benchmark solution described herein may be 
used for validation after porting to other platforms and dif- 
fering compilers. 35 

A similar analytic benchmark has been developed for the 
1977 Solar minimum galactic cosmic ray spectrum. This 
benchmark demonstrates that the propagator ignoring sec- 
ondary particle production and fragmentation are a fraction of 
percent of the corresponding analytic solution with main 40 
errors near the boundaries of the energy grid, as shown in 
FIGS, lla-b, and most values are correct to a small fraction of 
1 percent. The dose and dose equivalent of the analytic bench- 
mark solution and numerical benchmark solution differ by 
less than 0.15 percent. 45 

Benchmarking can be important in both evaluation of code 
accuracy as well as a provision of test cases for code verifi- 
cation after porting to other platforms and/or compilers. 

FIG. 12 is a plot 1200 illustrating a comparison of the 
results derived from the BRYTRN (version 3) 1260, the 1 995 50 
HZETRN 1262 (including ten years of drift), and the present 
method (improved numerical procedures as developed 
above). The plots 1200 shown in FIG. 12 demonstrate the 
differences in dose equivalent 1210 (ICRP 60 ) shielded at 55 
different depths of water 1220 from the Webber spectrum by 
20 g/cm 2 of iron between the different computations methods 
1260, 1262 and 1264. 

There are many reasons for the differences, including cor- 
ruption of a nuclear reaction routine for light ions and a 50 
nuclear fragmentation database, in addition to development 
of improved numerical procedures. Appropriate modifica- 
tions as discussed above have been made resulting in the 
present method having corrected nuclear routines and data- 
base. A benchmark was used based on the high-energy trans- 65 
port code (FIETC) result using the Webber spectrum of 30-cm 
slab of water shielded by 20 g/cm 2 iron (or aluminum), as 
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shown in FIG. 10 for the present method in comparison with 
dose and dose equivalent (ICRP 26 quality factor) according 
to HETC. 

The dose and dose equivalent in water are given in Table 2 
for 20 g/cm 2 shields of aluminum and iron herein below. 
Table 2 shows the dose (cGy) and dose equivalent (cSv) in a 
30 cm water slab protected by aluminum or iron shield from 
the Webber solar particle event spectrum. 


TABLE 2 


Water 
Depth, cm 

X 

Aluminum Shield 
Thickness of 20 g/cm 2 

Iron Shield 

Thickness of 20 g/cm 2 

D(x), cGy* 

H(x), cSv** 

D(x), cGy* 

H(x), cSv** 

0 

7.09 (6.83) 

11.86 (11.56) 

9.18 (8.84) 

15.39 

(15.12) 

5 

3.86 (3.75) 

6.06 (5.99) 

4.68 (4.54) 

7.32 

(7.26) 

10 

2.36 (2.28) 

3.84 (3.75) 

2.77 (2.68) 

4.45 

(4.37) 

15 

1.53 (1.48) 

2.53 (2.61) 

1.77 (1.71) 

2.95 

( 2 . 86 ) 

20 

1.04 (1.00) 

1.85 (1.79) 

1.18 (1.14) 

2.07 

(1.99) 

25 

0.74 (0.71) 

1.40 (1.32) 

0.83 (0.78) 

1.52 

(1.45) 

30 

0.54(0.51) 

1.08 ( 1 . 02 ) 

0.60 (0.57) 

1.16 

(1.09) 


*values in parentheses are expected for TLD100 
**values in parentheses are for ICRP 26 quality factors 


Values for the 1 977 Solar minimum GCR spectrum for the 
aluminum or iron shielded water are shown in Table 3. In 
Table 3, annual dose (cGy) and dose equivalent (cSv) in a 30 
cm water slab protected by aluminum or iron shield from the 
1977 Solar Minimum GCR spectrum. 


TABLE 3 


Water 
Depth, cm 

X 

Aluminum Shield 
Thickness of 20 g/cm 2 

Iron Shield 

Thickness of 20 g/cm 2 

D(x), cGy* 

H(x), cSv** 

D(x), cGy* 

H(x), cSv** 

0 

20.9 (18.9) 

76.0 ( 66 . 8 ) 

22.0 (19.7) 

85.5 (75.7) 

5 

19.0 (17.5) 

58.2 (51.7) 

19.4 (17.8) 

64.9 (57.5) 

10 

18.3 (17.0) 

51.2 (45.8) 

18.6 (17.3) 

55.8 (49.8) 

15 

17.7 (16.6) 

46.5 (41.9) 

18.1 (16.8) 

49.9 (44.7) 

20 

17.3 (16.2) 

43.3 (41.8) 

17.6 (16.4) 

45.9 (41.3) 

25 

16.9 (15.9) 

41.1 (37.2) 

17.2 (16.1) 

43.1 (39.9) 

30 

16.5 (15.5) 

39.4 (35.7) 

16.8 (15.8) 

41.0(37.1) 


*values in parentheses are expected for TLD100 
**valuses in parentheses are for ICRP 26 quality factors 


In Tables 2 and 3, values for dose, expected TLD100 
response, and dose equivalent with ICRP 26 and ICRP 60 
quality factors are given. 

Additional benchmarks are provided for the two shield 
configurations described above (20 g/cm 2 of iron or alumi- 
num shielding water) from the Monte Carlo Codes PHITS, 
general-purpose particle and heavy ion transport Monte Carlo 
code developed by the Japan Atomic Energy Agency (JAERI/ 
JAEA), and MULASSIS, a Geant4-based multilayered 
shielding simulation tool, developed by the European Space 
Agency (ESA). The Monte Carlo results for the Webber spec- 
trum shown in Table 4 are compared with data from the 
present method reproduced in Table 2. More particularly, 
Table 4 shows the dose (cGy) and dose equivalent (cSv) in a 
30-cm water slab protected by aluminum or iron shield from 
the Webber solar particle event spectrum evaluated using 
recent Monte Carlo codes PHITS and MULASSIS (in paren- 
theses). 
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TABLE 4 


20 


Water Aluminum Shield Iron Shield 

Depth (cm) Thickness of 20 g/cm 2 Thickness of 20 g/cm 2 


x D(x), cGy* H(x), cSv* D(x), cGy* H(x), cSv* 


0 

7.09 

(6.82 

± 1.3%) 

10.9 (10.67 h 

t 3.3%) 

9.21 

(8.95 h 

t 1.2%) 

14.6 

(14.12 : 

± 2.8%) 

5 

3.90 

(3.76 

± 1.8%) 

5.95 

(5.62 h 

t 4.8%) 

4.74 

(4.54 h 

t 1.5%) 

7.16 

(6.55 ± 

3.2%) 

10 

2.37 

(2.27 

± 2.2%) 

3.70 

(3.48 h 

t 7.2%) 

2.79 

(2.72 h 

t 2.0%) 

4.26 

(4.14 ± 

6.5%) 

15 

1.53 

(1.48 

± 2.8%) 

2.44 

(2.14 h 

t 6.3%) 

1.76 

(1.73 h 

t 2.5%) 

2.74 

(2.56 ± 

6.8%) 

20 

1.03 

(1.02 

±3.4%) 

1.70 

(1.62 d 

t 8.3%) 

1.17 

(1.15 H 

t 3.2%) 

1.87 

(1.80 ± 

8.9%) 

25 

.717 

(0.72 

± 4.3%) 

1.21 

(1.05 i 

t 7.0%) 

0.806 

(0.85 h 

t 3.8%) 

1.32 

(1.33 ± 

14.5%) 

30 

.511 

(0.51 

±5.3%) 

.843 

(0.87 h 

t 18.3%) 

0.565 

(0.60 h 

t 4.8%) 

0.902 

(0.94 ± 

9.9%) 


PHITS results for the 1 977 Solar Minimum GCR spectrum 
are given in Table 5. More particularly, Table 5 shows the 
annual dose (cGy) and dose equivalent (cSv) in a 30 cm water 
slab protected by aluminum or iron shield from the 1977 Solar 
Minimum GCR spectrum evaluated using the recent Monte 
Carlo codes. 


TABLE 5 


Water 
Depth, cm 

X 

Aluminum Shield 
Thickness of 20 g/cm 2 

Iron Shield 

Thickness of 20 g/cm 2 

D(x), cGy 

H(x), cSv 

D(x), cGy 

H(x), cSv 

0 

23.1 

69.9 

24.6 

83.9 

5 

22.0 

56.3 

22.5 

63.2 

10 

21.6 

49.2 

21.8 

53.3 

15 

21.2 

44.6 

21.3 

47.2 

20 

20.8 

41.1 

21.0 

43.1 

25 

20.3 

37.8 

20.4 

39.1 

30 

18.6 

32.6 

18.7 

33.5 


As can be seen, there are differences between deterministic 
and Monte Carlo approaches, which tend to grow near the exit 
of the water column and may be caused by neutron (and lesser 
proton) leakage on the back surface that is not present in the 
present method. There are other differences, especially for 
1 977 Solar Minimum GCR penetration problem, on the order 
of ten to twenty percent in dose and dose equivalent, but not 
exceeding operational requirements of ±30 percent. 

The present invention advances Green’s function methods 
to produce a method that is capable of being validated using 
high-energy ion beams, treats the off-axis scattering in the 
propagation of the light-ion/neutron propagator, uses march- 
ing procedures for forward produced components of the inter- 
actions, and evaluates the production source terms with broad 
angles with more appropriate angle dependent propagation 
techniques. Further, it provides a generalized method for 
three nonhomogeneous material regions that uses propaga- 
tors with higher-order local truncation errors. This can be 
readily recognized by comparing equation (41) as used in 
2005 HZETRN with equation (42) as used in 1 995 HZETRN, 
which allows improved control of error propagation in the 
basic marching procedures (see FIG. 12 , comparing line 1264 
with line 1260). The process for converting to dose and dose 
equivalent uses improved numerical procedures based on a 
ten point Gauss -Legendre formulation, which was not avail- 
able in 1995 HZETRN. The nuclear physics model for the 
absorption cross section calculations has also been revised 
from 1995 HZETRN. Moreover, analytical benchmarks are 
included for code verification and in Table 1 as a portable test. 
A benchmark with an early version of the Oak Ridge National 
Laboratory HETC Monte Carlo code is provided in the 
present method according to FIGS. li)a-b. Also, a benchmark 


using the present method is given in Tables 2 and 3 . Tables 4 
and 5 contain new Monte Carlo benchmarks for evaluation of 
Tables 2 and 3. 

FIG. 13 is a flow chart 1300 of an embodiment of the 
present invention. The main program and each subroutine or 
function module begins with a brief description of its pur- 
pose. The complete method 1300 consists of a HZETRN core, 
subroutines, and function modules. The method 1300 trans- 
ports galactic cosmic ray (GCR) particles in free space (geo- 
magnetic cutoffs are ignored) through a given thickness of the 
aluminum shield followed by a given depth of water. The 
HZETRN computation method 1300 includes an interface for 
providing input options 1310 . An environmental model data- 
base is provided as an input. The array dimensions for the 
energy grid points and isotope fragment numbers are also 
entered along with the year in the solar cycle that is to be used. 
Finally, the depth in the aluminum shield where dosimetric 
quantities are to be calculated is provided as an input. 

Data is provided to support the atomic and nuclear inter- 
actions 1320 . For the atomic interactions, the energy, range, 
and stopping-power database for water and aluminum are 
entered. For the nuclear interactions, the absorption and frag- 
mentation cross-section database for water and aluminum are 
entered. The step size for the numerical-analytical propaga- 
tion algorithm 1330 may be entered. Dosimetric quantities 
subroutine 1340 accepts quality factor specifications and 
alternate risk estimate approach specifications. The Dosim- 
etric quantities subroutine 1340 then calculates the dose and 
dose equivalent, which is the product of the input quality 
factor, Q, and the dose at a given point in human tissue. The 
output options 1350 include specifying the fluxes, doses, an 
alternate risk estimate and linear energy transfer (LET) spec- 
tra. The output of the present method 1300 may be phased in 
to complex geometry models for designing spacecraft radia- 
tion shields based on the output. 

FIG. 14 illustrates a system 1400 according to an embodi- 
ment of the present invention. Embodiments of the present 
invention may take the form of an entirely hardware embodi- 
ment, an entirely software embodiment or an embodiment 
containing both hardware and software elements. In a pre- 
ferred embodiment, the invention is implemented in software, 
which includes but is not limited to firmware, resident soft- 
ware, microcode, etc. Furthermore, embodiments of the 
present invention may take the form of a computer program 
product 1490 accessible from a computer-usable or com- 
puter-readable medium 1468 providing program code for use 
by or in connection with a computer or any instruction execu- 
tion system. 

For the purposes of this description, a computer-usable or 
computer readable medium 1468 can be any apparatus that 
can contain, store, communicate, propagate, or transport the 
program for use by or in connection with the instruction 
execution system, apparatus, or device. The medium 1468 
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may be an electronic, magnetic, optical, electromagnetic, 
infrared, or semiconductor system (or apparatus or device) or 
a propagation medium. Examples of a computer-readable 
medium include a semiconductor or solid state memory, mag- 
netic tape, a removable computer diskette, a random access 
memory (RAM), a read-only memory (ROM), a rigid mag- 
netic disk and an optical disk. Current examples of optical 
disks include compact disk-read only memory (CD-ROM), 
compact disk-read/ write (CD-R/W) and DVD. 

A system suitable for storing and/or executing program 
code will include at least one processor 1496 coupled directly 
or indirectly to memory elements 1492 through a system bus 
1420 . The memory elements 1492 can include local memory 
employed during actual execution of the program code, bulk 
storage, and cache memories which provide temporary stor- 
age of at least some program code in order to reduce the 
number of times code must be retrieved from bulk storage 
during execution. 

Input/output or I/O devices 1430 (including but not limited 
to keyboards, displays, pointing devices, etc.) can be coupled 
to the system either directly to the system or through inter- 
vening I/O controllers. 

Network adapters 1450 may also be coupled to the system 
to enable the system to become coupled to other data process- 
ing systems 1452 , remote printers 1454 or storage devices 
1456 through intervening private or public networks 1460 . 
Modems, cable modem and Ethernet cards are just a few of 
the currently available types of network adapters. 

Accordingly, the computer program 1490 comprise 
instructions which, when read and executed by the system 
1400 of FIG. 14 , causes the system 1400 to perform the steps 
necessary to execute the steps or elements of the present 
invention. For example, one embodiment of the system 1400 
calculates high-energy neutron/ion transport to a target of 
interest by performing operations that include storing data 
defining boundaries for a calculation of a high-energy neu- 
tron/ion transport to a target of interest; calculating the high- 
energy neutron/ion transport to the target of interest using 
numerical procedures selected to reduce local truncation 
error by including higher order terms and to allow absolute 
control of propagated error by ensuring truncation error is 
third order in step size, and using scaling procedures for flux 
coupling terms modified to improve computed results by 
adding a scaling factor to terms describing production of 
j -particles from collisions of k-particles; and providing the 
calculated high-energy neutron/ion transport to modeling 
modules to control an effective radiation dose at the target of 
interest. 

The foregoing description of the embodiment of the inven- 
tion has been presented for the purposes of illustration and 
description. It is not intended to be exhaustive or to limit the 
invention to the precise form disclosed. Many modifications 
and variations are possible in light of the above teaching. It is 
intended that the scope of the invention be limited not with 
this detailed description, but rather by the claims appended 
hereto. 

What is claimed as new and desired to be secured by 
Letters Patent of the United States is: 

1. A computer-implemented method for calculating a 
transport of a high-energy neutron/ion transport flux to a 
target of interest within a shielded region, comprising: 

defining boundaries for the transport of the high-energy 
neutron/ion transport flux to the target of interest within 
the shielded region; 
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receiving, as input, at least one of shielding dimensions, 
identification of shielding materials, high-energy neu- 
tron/ion flux at the boundaries, and a spatial location for 
the target of interest; 

5 calculating the transport of the high-energy neutron/ion 
transport flux to the target of interest via the equation 

%(x+h, r)=exp r+Vjh)+J. k (vj/v k )o Jk (r+Vjh/ 

2)^ k [x, r+(vj+v k )h/2]x$ 0 h exp{-Oj(r+Vjh/2)x '-a k 
[ r+{vj+v k )h/2){h-x')]})dx'=Qxv [-X>j(r,h)]%(x, r+ 

10 Vjh)+2. k (Vj/v k )Oj k (r +vfi.l 2)ty k [x,r+ (v frV k )h/2 ] x 

[ ex P { “ °/ ( r +Vjh/2)h} - exp { - o k [ r+(,Vj+v k )h/2]h }]/ 

{° k [ r + ( v ) + v*)/z/2)]-a / -(r +Vjk/2)}+0(h 2 ) 

wherein i^(x+h,r) andxp^x+l^r) are scaled fluxes for j -par- 
ticles and k-particles at an end of a subinterval compu- 
15 tational point, ^ and t, k are high-energy neutron/ion 
fluxes at the boundaries, x and h are spatial coordinates, 
r is a single residual range coordinate, v 7 and \ k are 
scaling factors associated with j -particles and k-par- 
ticles, respectively, and cr, and o k are cross-sections for 
20 the j -particles and k-particles, respectively; 

wherein, when calculating the high-energy neutron/ion 
transport flux to the target of interest, propagated error 
for values calculated by the computer implemented 
numerical method is controlled by controlling trunca- 
25 tion error as a third order in step size; 

wherein, when calculating the high-energy neutron/ion 
transport flux to the target of interest, the scaling factors 
are added to adjust for behavior associated with produc- 
tion of j-particles from collisions of k-particles; and 
30 wherein, when calculating the high-energy neutron/ion 
transport flux to the target of interest, the single residual 
range coordinate is introduced for all neutrons/ions in 
the computer implemented numerical method. 

2. The method of claim 1, wherein the scaling factor is 
35 defined by a ratio vjv k , wherein v i . i s Z^/Aj, v k is Z 7c 2 / A k , A is 

mass number and Z is charge number. 

3. The method of claim 1, wherein the calculating high- 
energy neutron/ion transport flux to the target of interest 
further comprises calculating high-energy neutron/ion trans- 

40 port flux through at least one shield material. 

4. The method of claim 1, wherein the calculating high- 
energy neutron/ion transport flux to the target of interest 
further comprises calculating high-energy neutron/ion trans- 
port flux to a selected tissue. 

45 5. The method of claim 1, wherein the calculating high- 

energy neutron/ion transport flux to the target of interest 
further comprises using a uniform grid distributed over two 
sub-domains to provide greater accuracy with less grid points 
than required by the Hilly uniform grid. 

50 6. The method of claim 5, wherein calculating the high- 

energy neutron/ion transport flux to the target of interest 
further comprises implementing a three-point Simpson’ s rule 
to reduce integration errors, when evaluating a number of 
j-particles resulting from collisions of k-particles, by using 
55 midpoint values of the improved interpolation with the uni- 
form grid distributed over two sub -domains. 

7 . The method of claim 1 , wherein the calculating high- 
energy neutron/ion transport flux to the target of interest 
further comprises adjusting a number of grid points to accom- 

60 modate simulation of geomagnetic cutoff effects while main- 
taining high numerical accuracy. 

8. The method of claim 1, wherein the calculating high- 
energy neutron/ion transport flux to the target of interest 
further comprises verifying accuracy of light-ion/neutron 

65 cross section routines. 

9 . The method of claim 1 , wherein calculating dosimetric 
quantities from the high-energy neutron/ion transport flux to 
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the target of interest further comprises implementing a ten- 
point Gauss-Legendre quadrature to improve correlation of 
the effective radiation dose to analytic results. 

10. The method of claim 1 further comprising validating 
the calculated high-energy neutron/ion transport flux using 5 
measured dosimetry and dynamic anisotropic environmental 
models. 

11. The method of claim 1, further comprising: 

calculating a dose from the flux of the high energy neutron/ 

ion transport to the target of interest. 10 

12. The method of claim 1, wherein the scaling factor 
corrects for light ion propagation associated with the produc- 
tion of j -particles from k-particles. 

13 . The method of claim 12, wherein the light ion particles 

comprise at least one of hydrogen or helium isotopes. 15 

14. The method of claim 1, wherein the single residual 
range coordinate comprises mapping, at low energies, for the 
high-energy neutron/ion transport flux to the target of inter- 
est. 

15. The method of claim 1, wherein calculating the trans- 20 
port of the high-energy neutron/ion transport flux to the target 

of interest is accomplished in steps from the boundaries to the 
target of interest. 

16. The method of claim 4, wherein the tissue is a tumor. 

17. The method of claim 10, wherein validating the calcu- 25 
lated high-energy neutron/ion transport flux using measured 
dosimetry and dynamic anisotropic environmental models 
occurs with respect to a predetermined vehicle design. 

18. A device configured to calculate a transport of a high- 
energy neutron/ion transport flux to a target of interest within 30 
a shielded region, comprising: 

memory for storing data defining boundaries for the trans- 
port of the high-energy neutron/ion transport flux to the 
target of interest within the shielded region; 

an input device for receiving, as input, at least one of 
shielding dimensions, identification of shielding mate- 
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rials, high-energy neutron/ion flux at the boundaries, 
and a spatial location for the target of interest; and 
a processor, coupled to the memory, for 

calculating the transport of the high-energy neutron/ion 
transport flux to the target of interest interest via the 
equation 

%(x+h, r)=exp [-X>j(r,h)]%(x, r+Vjh)+J. k {vj/v k )Gj k (r+Vjh/ 

2)%/x, r+(vj+v k )k/2]xf 0 h exp{-Oj(r+Vjh/2)x '-o k 
[ ^+(v / .+v^)/ 2 / 2 ) (/ 2 -X 1 )] })^'=exp [-tjinh^jix, r+ 
v,h)+lljXv/v k )o ik (r+v i hA 2)ip Jx, r+(v:+VjA)hf 2] x 
l ex P {~ a j( r +Vjk/2)k} - exp { - o k [ r+(vj+v k )h/2]h }]/ 

{°k[ r + ( v a ,-+v*)/2/2)]-CT / (r +Vjh/2)}+0(h 2 ) 

wherein i^.(x+h,r) andxp^(x+h,r) are scaled fluxes for j -par- 
ticles and k-particles at an end of a subinterval compu- 
tational point, ty and t, k are high-energy neutron/ion 
fluxes at the boundaries, x and h are spatial coordinates, 
r is a single residual range coordinate, v y . and \ k are 
scaling factors associated with j -particles and k-par- 
ticles, respectively, and and o k are cross-sections for 
the j -particles and k-particles, respectively, 
wherein, when calculating the high-energy neutron/ion 
transport flux to the target of interest, propagated error 
for values calculated by the corn cuter implemented 
numerical method is controlled by controlling trunca- 
tion error as a third order in step size, 
wherein, when calculating the high-energy neutron/ion 
transport flux to the target of interest, the scaling factor 
are added to adjust for behavior associated with produc- 
tion of j-particles from collisions of k-particles, and 
wherein, when calculating the high-energy neutron/ion 
transport flux to the target of interest, the single 
residual range coordinate is introduced for all neu- 
trons/ions in the computer implemented numerical 
method. 



